clear all
*cap log close
set more off
set seed 413

// other pars -------------
local boot = 250
local Y    = "cite_ny1"
// ------------------------

// construct bootstrap weights -------
use "${temp}/data_extrap", clear 
keep officerid 
duplicates drop 
gen bwt_0 = 1 
forval j = 1/`boot' {
	gen temp1 = rgamma(1,1)
	egen temp2 = sum(temp1)
	gen bwt_`j' = temp1/temp2
	drop temp1 temp2
}
qui compress
save "${temp}/bootwt", replace 
// -------------------------------------


// initialize bootstrap loop -----------
forval j = 0/`boot' {

	* setup data ------
	use "${temp}/data_extrap", clear 
	gen Y = `Y' 
	gen D = harsh
	merge m:1 officerid using "${temp}/bootwt", keep(1 3) nogen 
	gen bwt = bwt_`j'
	drop bwt_*

	* estimation via optimal bins -------
	binsreg Y Z [aw=bwt], absorb(totfe) nbins() savedata(bs) replace 
	preserve 
	use bs, clear 
	qui summ dots_fit if dots_binid==1 
	local y0 = `r(mean)'
	qui summ dots_binid 
	qui summ dots_fit if dots_binid==`r(max)'
	local y1 = `r(mean)'
	restore 

	* if first iteration, store bins for next file....
	if `j' == 0 {
		preserve 
		use bs, clear 
		save "${temp}/optimal_binscatter", replace
		restore
	}
	rm bs.dta 

	* get other parameters ----------------
	qui summ Y [aw=bwt]
	local y = `r(mean)'
	qui summ D [aw=bwt]
	local p = `r(mean)'
	reghdfe Y D [aw=bwt]
	local y0d0 = _b[_cons]
	local y1d1 = _b[_cons]+_b[D]

	* store parameters ------------------
	clear 
	set obs 1 
	foreach par in y p y0d0 y1d1 y0 y1 {
		gen `par' = ``par'' 
	}

	* add other parameters -------------
	gen ate = y1-y0 
	gen att = (y-y0)/p 
	gen atu = (y1-y)/(1-p)
	gen y1d0 = (1/(1-p))*y1-(p/(1-p))*y1d1
	gen y0d1 = (1/p)*y0-((1-p)/p)*y0d0
	gen y0d1my0d0 = y0d1-y0d0
	gen attmatu = att-atu 

	* store bootstrap output ------------
	gen iter = `j'
	if `j'==0 {
		save "${temp}/nonparboot", replace
	}
	if `j'>=1 {
		qui append using "${temp}/nonparboot" 
		save "${temp}/nonparboot", replace 
	}

	* progress report ---------------
	di "Iteration `j' of `boot' completed"

}


// compile bootstrap output ------------------
// annoying stuff to format output same as R output
use "${temp}/nonparboot", clear 
keep if iter > 0 
collapse (sd) y-attmatu
tempfile se
save    `se'

use "${temp}/nonparboot", clear 
keep if iter == 0 
drop iter 
tempfile est 
save    `est'

ds
local varlist `r(varlist)'
local parcount: word count `varlist'

clear
set obs `parcount' 
gen par = ""
gen est = . 
gen se  = .
save "${est}/optimal_bins", replace  

local k=1 
foreach v in `varlist' {

	di "`v'"
	use `est', clear 
	summ `v'
	local est_copy = `r(mean)'
	use `se', clear 
	summ `v'
	local se_copy  = `r(mean)'

	use "${est}/optimal_bins", clear 
	replace par = "`v'" if _n == `k'
	replace est = `est_copy' if _n == `k'
	replace se  = `se_copy'  if _n == `k'
	save "${est}/optimal_bins", replace
	local ++k

}

gen ub = est+1.96*se 
gen lb = est-1.96*se
save "${est}/optimal_bins", replace
rm   "${temp}/nonparboot.dta"
rm   "${temp}/bootwt.dta"











